MacFormat España 14/MacFormat n. 14 (Spain)/MacFormat 14.bin/Shareware Internet/Educación/Ear Workout 2.1/source code/get_period.cp

Text File | 1996-01-09 | 15.2 KB | 431 lines

    int
get_period(
    unsigned char *data,    //-- your data
    int n,            //-- number of data points
    int lo_p,        //-- shortest possible period
    int hi_p,        //-- longest possible period
    short *grid,        //-- storage you allocate for me (call me with n=0
                //    to find out how many bytes)
    int grid_space,        //-- number of bytes you allocated for me
    int *naverage,        //-- number of periods averaged together to get this
                //    result (output)
    short **nhit_ptr,    //--pointer to the histogram; 256 bins per octave;
                //    bin number=256*log2(p/lo_p)
    int *nhist,        //--tells them how many bins in histogram,
    int **int_info_handle,    //--Passes back a pointer to some handy data:
                //    *int_info[0]=mean abs. val of deviation from mean
                //    *int_info[1]=mean abs. val of derivative
                //    (Means are only based on a subset of the
                //    data.)
                //    *int_info[2]=xshift
                //    *int_info[3]=yshift
                //    *int_info[4]=tstep
    double **double_info_handle
                //
    );

/*
Function get_period returns period of the waveform in units of the
    sampling period divided by 256.
Returns 0 on error.
If you call it with n=0, it returns the number of bytes you should
allocate for it storage.

This is intended to work for a slightly non-periodic waveform, i.e. one
whose fourier spectrum consists of peaks at frequencies which are very nearly
multiples of a fundamental frequency.

The concept:
One could just watch the waveform y(t) pass through a certain point (say y=0)
and then watch for the next time it passes throught that point.  The difference
in time would be the period.  Two problems:
(1) Periodic waveforms may pass through the same y value several times per period.
(2) For a typical sampling frequency, one period may be only 10-100 samples, so
    we'd like to average several measurements together.

Solution to (1):
A truly periodic waveform would come back to the same y', y'', y''', ...,
not just the same y.  We can make a grid in the "phase space" that has
y on one axis and y' on the other axis, and watch for when the wavform comes
back to the same cell in phase space.  As a matter of fact, a curve in a bounded
2-d region pretty much has to cross itself, so we're almost guaranteed to hit
cells repeatedly.  (We could generalize this to 3-d, but the property of guaranteed
crossing is lost, and also the 2nd derivative is very granular with an 8-bit
adc.)
Looking at phase space plots of almost periodic functions shows that
a typical plot is sort of a cardiod shape, possibly with a little loop in it
that is not a full period.  However, the self-crossing created by the little
loop does not usually have the same 2nd derivative -- the curve typically executes
that crossing at a near-right angle.  The self-crossing at one full period, however,
typically happens at a shallow angle.  This means that the correct full-period
crossing generates a lot of return visits to cells in phase space, while the
incorrect crossing generated by the little loops generates far fewer return
visits.  So a good method is to look for lots of crossings that have about the
same period.

Solution to (2):
Every time we get a crossing, take the log base 2 of its period (times 256).
Make a histogram of this quantity.  The peaks of the histogram are the
periods we've seen frequently.  For each histogram bin, also store the sum of
the periods, so we can divide later by the number of hits to get an average period.
Also include neighboring histogram bins in the average to avoid funny effects
associated with bin boundaries.

Fine points:

How to decide on the size of the cells in phase space?
Take a sample from the data and take the average absolute value.  Do a bit shift
on the y and y' values so that the resulting range approximately fills half the
size of the grid.  If we make our cells in phase space too small, we may not
get any crossings -- if this happens, increase the bit shift by one and try again.

We may have several data points in a row in the same cell of phase space.  Only
use the first one to fall into a given cell.

Actual samples have a tendency to drift across the phase space plot -- I subtract
out a running average from each data point.

Current status:
Tested with singing, whistling and guitar.  Works well for signals of sufficient
amplitude.  For weak signals, the granularity is too poor for the phase-space
concept to work, and the method tends to find a broad, strong false histogram peak
at low periods and only a much weaker peak at the correct, longer period.  For
most of these weak signals, though, I can still hear a note with a definite pitch
when I play it back.  This points up the advantages of the frequency domain:
you get a sum over many periods, which improves the signal-to-noise.
Seems to work quite well with whistling, which is probably more sinusoidal,
and not as well with singing.  With singing, seems to get confused by higher
harmonics.  Do some integration?

Idea:
Normally you want some _low-pass_ filtering before you do your pitch detection,
but by using y and y' for my phase-space plots, I'm essentially doing some
_high-pass_ filtering (differentiation to get y').  How about constructing the
phase-space plot from y and the _integral_ of y?  Would have to blank it
out at start to avoid transients.

Possible improvement: First do the phase-space method,
then fft with severe resampling to make it fast enough.
The phase-space method is immune to
aliasing (other than that of the hardware sampling freq), but prone
to the problem described above for weak signals.  The fft is more
likely to get a good peak for a weak signal, but has problem with aliasing.
If the fft and phase-space method agree on some fairly strong peak, that's
probably it.

Problem:  Tries to extract a period, even if the sound is just noise (e.g.
the consonant "sh").

Notes from Musical Applications of Microprocessors, Hal Chamberlin, Hayden, 1980:
  A big challenge: sounds with _quick decays_, because they're non-periodic.
  Suggested test suite:
      a) sine wave
      b) sine-like wave with very flat crests
      c) missing fundamental
      d) strong harmonic
      e) wave in which a given y and y' are repeated every half-period
      f) speech
      g) rapid tremolo
      h) white noise (should be correctly flagged by routine)
  Time-domain methods:
      low-pass filtering to emphasize fundamental, but:
          repeated integration emphasizes low-f transients, and
              may blank out routine for a long time
          don't know a priori where fundamental is
          might not have fundamental
      my idea: would distortion restore a missing fundamental?
      typical method:
          Put signal through diode to select only pos peaks.
          Charge up capacitor on peak of waveform.  Attenuate 1-5%.
          Define period in terms of crossing of original and
          altered signal.
      variation:
          Do the above once with positive side of signal, once
          with neg side.  If 2 estimates disagree, choose the
          lower period.
      Gold & Rabiner:
          Not really clearly described.  Six-way majority logic,
          involving peak-peak, valley-valley, peak-valley,...
          Decision tree using these data combined with result
          of previous two determinations.
  Autocorrelation:
      Can use true autocorrelation, f(x) * f(x+k), or things
      like |f(x)-f(x+k)|.  Look for max or min of autocorrel w.r.t. k.
      Weaknesses: can give peaks at other harmonics; pitch shifts
       if waveform changes.
  Cepstrum techniques:
      Def of cepstrum:
          Do fft to get power spectrum, take log on y axis.
          Take another forward (?) fft.
          This gives a cepstrum plot.  X axis has dimensions of time,
          but it's called quefrency.
      Idea: power spectrum has a bunch of harmonics, evenly spaced.
          Fft on that gives spacing of harmonics, which is
          freq. of fundamental.
          Low-quefrency part of plot tells about envelope of resonances,
          high-quefrency tells about the harmonics that excited
          the resonances.
          Reason: by taking log of (excitation fn)*(resonance curve),
          you're making it into a _sum_, which fft can separate.
      My misgiving:  what if power spectrum is very small at some point...
          then log power spectrum has huge neg peak?????
      Problem: fails if only odd harmonics are present.
 */


unsigned char log2_256_lookup(unsigned char x);
int quick_log2_ratio(int x,int y);

#define XGRIDSHIFT 10
#define XGRID (1<<XGRIDSHIFT)
#define YGRID 128
#define TOTGRID (XGRID*YGRID)
#define ABSVAL(x) ((x)>=0 ? (x) : -(x))
#define MAXOCTAVES 20
            //...max of 20 octaves between lo_p and hi_p
#define PREC 8
    //...number of bits of precision of input data

#define DEBUG 0

#if DEBUG
#include <stdio.h>
#endif

    int
get_period(
    unsigned char *data,    //-- your data
    int n,            //-- number of data points
    int lo_p,        //-- shortest possible period
    int hi_p,        //-- longest possible period
    short *grid,        //-- storage you allocate for me (call me with n=0
                //    to find out how many bytes)
    int grid_space,        //-- number of bytes you allocated for me
    int *naverage,        //-- number of periods averaged together to get this
                //    result (output)
    short **nhit_ptr,    //--pointer to the histogram; 256 bins per octave;
                //    bin number=256*log2(p/lo_p)
    int *nhist,        //--tells them how many bins in histogram,
    int **int_info_handle,    //--Passes back a pointer to some handy data:
                //    *int_info[0]=mean abs. val of deviation from mean
                //    *int_info[1]=mean abs. val of derivative
                //    (Means are only based on a subset of the
                //    data.)
                //    *int_info[2]=xshift
                //    *int_info[3]=yshift
                //    *int_info[4]=tstep
    double **double_info_handle
                //
    )
    {
      static short nhit[MAXOCTAVES*256];
      static short sum_of_periods[MAXOCTAVES*256];
      static int int_info[20];
      static double double_info[20];
      int i,j,x,y,u,v,xnorm,ynorm,nnorm,xshift,yshift,q,p,old,log2p,any_hits,
          most_popular,log2p_below,log2p_above,nvotes,this_period,
          best_period,most_votes,old_cell,running_sum,running_sum_time,
          fake_running_sum_value,k,running_average,general_average,
          tstep;

      if (nhit_ptr != (short **) 0) *nhit_ptr = nhit;
      if (nhist != (int *) 0) *nhist = MAXOCTAVES*256;
      if (int_info_handle != (int **) 0) *int_info_handle = int_info;
      if (double_info_handle != (double **) 0) *double_info_handle = double_info;

      if (n==0) return TOTGRID*sizeof(short);

      if (grid_space<TOTGRID*sizeof(short)) return 0;

      //-- get a fairly random sample of
      //   typical magnitude of x & y, and determine DC offset
      xnorm = 0;
      ynorm = 0;
      nnorm = 0;
      general_average = 0;
      q = lo_p;
      while (n/q<10) q *= 2;
      for (i=1; i<n; i+=q) {
        general_average += data[i];
        ++nnorm;
      }
      general_average /= nnorm;
      for (i=1; i<n; i+=q) {
        u = data[i]; //--convert data to signed
        v = data[i-1];
        x = u-general_average;
        y = u-v;
        xnorm += ABSVAL(x);
        ynorm += ABSVAL(y);
      }

      xnorm /= nnorm;
      ynorm /= nnorm;
      int_info[0] = xnorm;
      int_info[1] = ynorm;

      xshift = 0;
      while (xnorm>>xshift > XGRID/4)
        ++xshift;

      yshift = 0;
      while (ynorm>>yshift > YGRID/4)
        ++yshift;
      if (yshift==0 && ynorm>0 && ynorm < YGRID/8)
        tstep = YGRID/(8*ynorm);
      else
        tstep = 1;
      if (tstep<1) tstep=1;
      if (tstep>lo_p/4) tstep=lo_p/4;

      int_info[2] = xshift;
      int_info[3] = yshift;
      int_info[4] = tstep;

      fake_running_sum_value = general_average;
      running_sum_time = lo_p/2;

      for (;;) {
          any_hits = 0;
          for (i=0; i<TOTGRID; i++) {
            grid[i] = -1;
          }
          for (i=0; i<MAXOCTAVES*256; i++) {
            nhit[i] = 0;
            sum_of_periods[i] = 0;
          }
          running_sum = fake_running_sum_value*running_sum_time;
          for (i=tstep; i<n; i++) {
            running_sum += data[i];
            k = i-running_sum_time;
            if (k>0)
              running_sum -= data[k];
            else
              running_sum -= fake_running_sum_value;
            running_average = running_sum/running_sum_time;
            u = data[i]; //--convert data to signed
            v = data[i-tstep];
            x = u-running_average;
            y = u-v;
            x = x>>xshift + XGRID/2;
            if (x<0) x=0;
            if (x>XGRID-1) x=XGRID-1;
            y = y>>yshift + YGRID/2;
            if (y<0) y=0;
            if (y>YGRID-1) y=YGRID-1;
            j = x + y<<XGRIDSHIFT;
            old = grid[j];
            grid[j] = i;
            if (old>=0 && j!=old_cell && i!=1 ) {//-- deja vu
              p = (i-old); //-- possible period
              if (p>=lo_p && p<=hi_p) {
                log2p = quick_log2_ratio((int) p,(int) lo_p);
                if (log2p>=0 && log2p<MAXOCTAVES*256) {
                  ++nhit[log2p];
                  sum_of_periods[log2p] += p;
                  any_hits = 1;
                }
              }
            }
            old_cell = j;
          }
          if (any_hits) {
            //-- find out which period is the most popular
            most_votes = 0;
            for (p=lo_p; p<=hi_p; p++) {
              log2p_below = quick_log2_ratio((int) (p-1),(int) lo_p);
              log2p_above = quick_log2_ratio((int) (p+1),(int) lo_p);
              nvotes = 0;
              this_period = 0;
              for (log2p=log2p_below; log2p<=log2p_above; log2p++) {
                if (log2p>=0 && log2p<MAXOCTAVES*256) {
                  nvotes += nhit[log2p];
                  this_period += sum_of_periods[log2p];
                }
              }
              if (nvotes>most_votes) {
                best_period = 256*this_period;
                best_period = best_period/nvotes;
                most_votes = nvotes;
              }
            }
            if (most_votes>0) {
              *naverage = most_votes;
              return best_period;
            }
          }
          if (xshift>=PREC-1 || yshift>=PREC-1) return 0; //-- we failed;
          xshift++;
          yshift++;
      }



    }

    //-- returns 256*(log base 2 of x/y)
    // returns garbage (result=32000 or -32000) if x<=0 or y<=0
    int
quick_log2_ratio(int x,int y)
    {
      int z,q;
      unsigned char qq;
      if (y<=0) return 32000;
      if (x<=0) return -32000;
      z = 0;
      while (x>=2*y) {
        y = y*2;
        z += 256;
      }
      while (x<y) {
        x = x*2;
        z -= 256;
      }
      q = 256*(x-y);
      q = q/y;
      if (q<=255)
        qq = q;
      else
        qq = 255;
      z += log2_256_lookup(qq);
      return z;
    }


/*
    Made the lookup table in the function log2_256_lookup() using
    the following code:
          {
            FILE *f;
            int i;
            double x;
            f = fopen("a.a","w");
            for (i=0; i<=255; i++) {
              x = i/256.;
              fprintf(f,"%d",(int) (log(1.+x)/log(2.)*256.));
              if (i<255) fprintf(f,",");
              if ((i+1)%20==0) fprintf(f,"\n");
            }
            fclose(f);
          }
*/
    unsigned char
log2_256_lookup(unsigned char x)
    {
      static unsigned char table[256] = {
0,1,2,4,5,7,8,9,11,12,14,15,16,18,19,21,22,23,25,26,
27,29,30,31,33,34,35,37,38,39,40,42,43,44,46,47,48,49,51,52,
53,54,56,57,58,59,61,62,63,64,65,67,68,69,70,71,73,74,75,76,
77,78,80,81,82,83,84,85,87,88,89,90,91,92,93,94,96,97,98,99,
100,101,102,103,104,105,106,108,109,110,111,112,113,114,115,116,117,118,119,120,
121,122,123,124,125,126,127,128,129,131,132,133,134,135,136,137,138,139,140,140,
141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,
161,162,162,163,164,165,166,167,168,169,170,171,172,173,173,174,175,176,177,178,
179,180,181,181,182,183,184,185,186,187,188,188,189,190,191,192,193,194,194,195,
196,197,198,199,200,200,201,202,203,204,205,205,206,207,208,209,209,210,211,212,
213,214,214,215,216,217,218,218,219,220,221,222,222,223,224,225,225,226,227,228,
229,229,230,231,232,232,233,234,235,235,236,237,238,239,239,240,241,242,242,243,
244,245,245,246,247,247,248,249,250,250,251,252,253,253,254,255};
      return table[x];
    }